
# 1.加载所需库------
if (!require("pacman")) install.packages("pacman")
pkgs <- c("dplyr", "ggplot2", "readxl","writexl", "tidyr", "lubridate", "janitor", "stringr", "purrr",
          "psych", "magrittr") #, "sf", "spdep", "tmap", "ggspatial"
pacman::p_load(pkgs, character.only = TRUE)

# 2.读取人口数据,排除兵团------
#读取中英文名字
nameCE <- read_excel("data/origin/发病率2004-2023年v1.xlsx", 
                     sheet = "中英文名") %>%
  dplyr::select(region,county_code,city_code,nameE,nameC = 地区)
nameCE$county_code %>% unique() %in% "654203" %>% sum()
nameCE$county_code %>% unique() %in% "654223" %>% sum()
nameCE$county_code <- recode(nameCE$county_code,"654203" = "654223")


nameCE_city <- nameCE %>% filter(region == "city") %>% 
  dplyr::select(city_code, nameE) %>% distinct()

nameCE_county <- nameCE %>% filter(region == "county") %>% 
  dplyr::select(county_code, nameE) %>% distinct()
nameCE$county_code %>% unique() %in% "654223" %>% sum()

# 生成省、市、区县的年、月的时空面板数据

# 生成时空面板数据
spacetime_panel <- nameCE %>%
  crossing(
    tibble(
      date = seq.Date(
        from = as.Date("2005-01-01"),
        to = as.Date("2023-12-01"),
        by = "month")) %>%
      mutate(
        year = year(date),
        month = month(date),
        date = format(date, "%Y/%m")))   # 标准化日期格式
      
# 结果验证
glimpse(spacetime_panel)

spacetime_panel$county_code %>% unique() %in% "654223" %>% sum()

# 读取人口数据
library(readxl)
pop <- read_excel("data/origin/发病率2004-2023年v1.xlsx", 
       sheet = "Report(基于六普七普人口数)", col_types = c("text", "text", "text", 
       "skip", "skip", "numeric", "numeric", "numeric", "numeric", "numeric", 
       "numeric", "numeric", "numeric","numeric", "numeric", "numeric", "numeric", 
       "numeric", "numeric","numeric", "numeric", "numeric", "numeric", "numeric", "skip")) %>%
  filter(region != "county兵团") 

pop_long <- pop %>% pivot_longer(cols = -c(地区,地区编码,region),names_to = "year") %>%
  mutate(
    地区编码 = case_when(
      region == "city" ~ str_sub(地区编码, 1, 4),  # 取前4位[1](@ref)
      region == "county" ~ str_sub(地区编码, 1, 6),  # 取前6位[1](@ref)
      TRUE ~ 地区编码  # 其他情况保持原值
    )
  )


pop_long$地区编码 %>% unique()

pop_long$地区编码 %>% unique() %in% "654223" %>% sum()#FALSE
pop_long$地区编码 %>% unique() %in% "654203" %>% sum()# TRUE

pop_long$地区编码 <- recode(pop_long$地区编码, `654203` = '654223') #shawan
pop_long$地区编码 %>% unique() %in% "654223"  %>% sum()#TRUE

pop_long$地区编码 %>% unique(); 96+14+1
pop_long$地区编码 %>% unique() %in% "654223" %>% sum()
pop_long$region %>% unique()
pop_long %>% dplyr::filter(region == "city");14*19
pop_long %>% dplyr::filter(region == "county") %>% unique();96*19



# 读取病例数据
cleaned_data9 <- read_excel("./data/processed/cleaned_data9.xlsx") %>%
  filter(year >= 2005)

cleaned_data9$city_code %>% unique() %>% length()
cleaned_data9$county_code %>% unique()%>% length();95
cleaned_data9$county_code %>% unique() %in% "654223" %>% sum()


# 3. 发病率------
#3.1逐年发病率------

# 1. 按全省分组的每年的病例数
province_level_cases <- cleaned_data9 %>%
  group_by(year) %>%
  summarise(count = n(), .groups = 'drop') %>%
  mutate(category = "province", subgroup = "650000") %>%
  dplyr::select(category,subgroup,year,count)

dim(province_level_cases)

# 2. 按区县（county_code）分组的每年的病例数
spacetime_panel_county <- spacetime_panel %>% filter(region == "county") %>% 
  dplyr::select(category = region, subgroup = county_code,year) %>% 
  distinct()
spacetime_panel_county %>% filter(subgroup == "654223")

county_level_cases0 <- cleaned_data9 %>%
  group_by(year, county_code) %>%
  summarise(count = n(), .groups = 'drop') %>%
  mutate(category = "county", subgroup = county_code) %>%
  dplyr::select(category, subgroup, year, count)
county_level_cases <- county_level_cases0 %>% right_join(spacetime_panel_county) %>% 
  mutate(count = coalesce(count, 0L))  # 0L表示整数型零 #mutate(count = ifelse(is.na(count), 0, count))

spacetime_panel_county %>% filter(subgroup == "654223")
county_level_cases %>% filter(subgroup == "654223")


# 3. 按城市（city_code）分组的每年的病例数
spacetime_panel_city <- spacetime_panel %>% filter(region == "city") %>% 
  dplyr::select(category = region, subgroup = city_code,year) %>% 
  distinct()
city_level_cases0 <- cleaned_data9 %>%
  group_by(year, city_code) %>%
  summarise(count = n(), .groups = 'drop') %>%
  mutate(category = "city", subgroup = city_code) %>%
  dplyr::select(category, subgroup, year, count)
city_level_cases <- city_level_cases0 %>% right_join(spacetime_panel_city) %>% 
  mutate(count = coalesce(count, 0L))
dim(city_level_cases)

# 合并所有结果,排除兵团的
case_year_long <- bind_rows(
  province_level_cases,
  county_level_cases,
  city_level_cases) %>% filter(category != "county兵团")
dim(case_year_long)

# 将数据转为宽格式以便展示
case_year <- case_year_long %>%
  pivot_wider(names_from = year, values_from = count)#, names_prefix = "Y"
dim(case_year)
dim(pop)
names(case_year);names(pop)
case_year %>% filter(category =="city")

# 保存人口和病例数据
openxlsx::write.xlsx(case_year, "./data/processed/case_year.xlsx")
openxlsx::write.xlsx(pop, "./data/processed/pop_year.xlsx")

# 整合人口和病例数据
pop_long
case_year_long

# 确保地区编码为前6位，并转换year为数值型
pop_long <- pop_long %>%
  mutate(地区编码 = substr(地区编码, 1, 6), year = as.double(year))

# 合并数据集，使用0填充NA值
rate_year <- pop_long %>%
  left_join(case_year_long, by = c("region" = "category", "地区编码" = "subgroup", "year" = "year")) %>%
  mutate(across(c(count), ~coalesce(.x, 0))) %>% 
  mutate(rate = count/value*100000)
  
# 查看结果
print(rate_year)
rate_year_province <- rate_year %>% filter(region == "province");19;dim(rate_year_province)
rate_year_city <- rate_year %>% filter(region == "city");14*19;dim(rate_year_city)
rate_year_county <- rate_year %>% filter(region == "county");96*19;dim(rate_year_county)
# rate_year_county %>% View()
# # 保存逐年发病率数据
# openxlsx::write.xlsx(rate_year, "./data/processed/rate_month.xlsx")

#3.2分阶段全省发病率------

periods <- list(
  # "All Years" = range(data$year),
  "2005-2009" = c(2005, 2009),
  "2010-2015" = c(2010, 2015),
  "2016-2019" = c(2016, 2019),
  "2020-2023" = c(2020, 2023)
)

# 计算各阶段统计量:累计发病率
result <- bind_rows(lapply(names(periods), function(p) {
  start <- periods[[p]][1]
  end <- periods[[p]][2]
  
  rate_year %>%
    filter(region == "province",        # 筛选省级数据[4](@ref)
           between(year, start, end)) %>% # 时间段筛选[2,4](@ref)
    summarise(
      时间段 = p,
      发病数 = sum(count, na.rm = TRUE),
      总人口数 = sum(value, na.rm = TRUE),
      发病率 = (发病数 / 总人口数) * 100000  # 按十万人口计算[4,5](@ref)
    )
}))

# 计算各阶段统计量：各阶段平均年发病率
result <- map_dfr(names(periods), function(p) {
  start <- periods[[p]][1]
  end <- periods[[p]][2]
  
  stage_data <- rate_year %>%
    filter(region == "province",        # 筛选省级数据
           between(year, start, end))  # 时间段筛选
  
  total_count <- sum(stage_data$count, na.rm = TRUE)
  total_pop <- sum(stage_data$value, na.rm = TRUE)
  n_years <- end - start + 1           # 计算实际年数
  
  data.frame(
    时间段 = p,
    年平均发病数 = round(total_count / n_years, 1),
    年平均发病率 = round((total_count / total_pop) * 100000, 1)
  )
})

# 结果格式化输出
result %>% 
  mutate(年平均发病率 = format(年平均发病率, nsmall = 1)) %>% 
  knitr::kable(align = "c")
# 结果展示
print(result)


#3.3分阶段各地州、区县发病率------

library(dplyr)
library(openxlsx)

periods0 <- list(
  "2005-2009" = 2005:2009,
  "2010-2015" = 2010:2015,
  "2016-2019" = 2016:2019,
  "2020-2023" = 2020:2023
)

# 添加阶段标签到原始数据
rate_year1 <- rate_year %>%
  mutate(
    period = case_when(
      year %in% periods0[[1]] ~ names(periods0)[1],
      year %in% periods0[[2]] ~ names(periods0)[2],
      year %in% periods0[[3]] ~ names(periods0)[3],
      year %in% periods0[[4]] ~ names(periods0)[4]
    )
  )

分阶段统计 <- rate_year1 %>%
  group_by(地区, 地区编码, region, period) %>%
  summarise(
    # 累计发病数
    累计发病数 = sum(count, na.rm = TRUE),
    # 年均发病数 = 累计发病数 / 阶段年数
    年均发病数 = 累计发病数 / length(unique(year)),
    # 年均人口 = 阶段内各年人口均值
    年均人口 = mean(value, na.rm = TRUE),
    # 年均发病率 = (年均发病数 / 年均人口) × 10^5
    年均发病率 = round((年均发病数 / 年均人口) * 1e5, 2)
  ) %>%
  ungroup()

openxlsx::write.xlsx(分阶段统计, "./output/发病率_分阶段.xlsx")
分阶段统计$region %>% unique()

#4.逐年逐月发病率------

# 1. 按全省分组的每年的病例数
spacetime_panel_province <- spacetime_panel %>% filter(region == "province") %>% 
  dplyr::select(category = region, subgroup = county_code,year, month) %>% 
  distinct()

province_level_cases0 <- cleaned_data9 %>%
  group_by(year,month) %>%
  summarise(count = n(), .groups = 'drop') %>%
  mutate(category = "province", subgroup = "650000") %>%
  dplyr::select(category,subgroup,year,month,count)
province_level_cases <- province_level_cases0 %>% right_join(spacetime_panel_province)%>% 
  mutate(count = coalesce(count, 0L))  ;1*19*12
dim(province_level_cases)

# 2. 按区县（county_code）分组的每年的病例数
spacetime_panel_county <- spacetime_panel %>% filter(region == "county") %>% 
  dplyr::select(category = region, subgroup = county_code,year, month) %>% 
  distinct()

county_level_cases0 <- cleaned_data9 %>%
  group_by(year, month,county_code) %>%
  summarise(count = n(), .groups = 'drop') %>%
  mutate(category = "county", subgroup = county_code) %>%
  dplyr::select(category, subgroup, year,month, count)
county_level_cases <- county_level_cases0 %>% right_join(spacetime_panel_county)%>% 
  mutate(count = coalesce(count, 0L))  
dim(county_level_cases);96*19*12

# 3. 按城市（city_code）分组的每年的病例数
spacetime_panel_city <- spacetime_panel %>% filter(region == "city") %>% 
  dplyr::select(category = region, subgroup = city_code,year, month) %>% 
  distinct()

city_level_cases0 <- cleaned_data9 %>%
  group_by(year, month,city_code) %>%
  summarise(count = n(), .groups = 'drop') %>%
  mutate(category = "city", subgroup = city_code) %>%
  dplyr::select(category, subgroup, year, month,count)
city_level_cases <- city_level_cases0 %>% right_join(spacetime_panel_city)%>% 
  mutate(count = coalesce(count, 0L))  

dim(city_level_cases); 14*12*19

# 合并所有结果,排除兵团的
case_month_long <- bind_rows(
  province_level_cases,
  county_level_cases,
  city_level_cases) %>% filter(category != "county兵团")
dim(case_month_long)

# 将数据转为宽格式以便展示
# case_month <- case_month_long %>%
#   pivot_wider(names_from = c(year,month), values_from = count)#, names_prefix = "Y"
# dim(case_month)
dim(pop)


# # 保存人口和病例数据
# openxlsx::write.xlsx(case_year, "./data/processed/case_year.xlsx")
# openxlsx::write.xlsx(pop, "./data/processed/pop_year.xlsx")

# 创建月份列
months <- tibble(month = 1:12)

# 使用crossing函数生成所有组合，并合并数据
pop_long_expanded <- pop_long %>%  crossing(months);(1+14+96)*19;dim(pop_long)



# 整合人口和病例数据
pop_long_expanded;dim(pop_long)[[1]]*12;2109*12;dim(pop_long_expanded)
case_month_long

# 确保地区编码为前6位，并转换year为数值型
pop_long_expanded <- pop_long_expanded %>%
  mutate(地区编码 = substr(地区编码, 1, 6), year = as.double(year))

# 合并数据集，使用0填充NA值
rate_month <- pop_long_expanded %>%
  left_join(case_month_long, 
            by = c("region" = "category", "地区编码" = "subgroup", "year", "month")) %>%
  mutate(across(c(count), ~coalesce(.x, 0))) %>% 
  mutate(rate = count/value*100000)

# 查看结果
print(rate_month)

rate_month_rate <- rate_month%>% filter(region == "province") %>%
  


# 保存逐月发病率数据
openxlsx::write.xlsx(rate_month, "./data/processed/rate_month.xlsx")


